clc
close all
clear all

D = 3
X = rand(D, D);

S = inverseObjective(X);
gS = inverseGradient(X);

for i = 1:D
    for j = 1:D
        X2 = X;
        X2(i, j) = X2(i, j) - 1e-6;
        X3 = X;
        X3(i, j) = X3(i, j) + 1e-6;
        gg(:, :, i, j) = 0.5*(inverseObjective(X3)-inverseObjective(X2))/1e-6;
    end
end

sum((gg(:)-gS(:)).^2)